In the past few decades, the study of binary compounds under high pressure has attracted significant attention from scientists and researchers thanks to their potential excellent properties, such as high-temperature superconductivity, extreme hardness, and antimetallization.^[1–8] Boride has been widely studied and applied in recent years. The cubic boron nitride (c-BN) shows excellent hardness, which is considered to be a promising superhard material in both technological and industrial applications.^[9] A superconducting (SC) critical temperature T_c of 39 K has been observed in MgB₂.^[10] Additionally, Fe₂B consisting of 99.97 wt% pure iron was found to provide excellent hardness and fracture toughness.^[11] Titanium diboride has been applied in many fields due to some interesting properties, such as high melting point, excellent electrical conductivity, and thermal conductivity.^[12]

Boron is a fascinating and complex element with only three valence electrons in the valence orbital, which is more likely to lose electrons and then become a common cation. In a room environment, pure boron and most of its compounds can be classified as insulators and semiconductors. Pressure can remarkably change the physical properties of materials and induce some unexpected and interesting physical phenomena. In the last few years, B has been widely investigated and was proposed to possess several fascinating structures, such as α-B₁₂ (space group R-3m), β-B₂₈ (space group Pnnm), and α-Ga (space group Cmca).^{[13, 14]} The phase of α-B₁₂ icosahedron is close to β-B₁₀₆ in terms of energy at ambient conditions, and becomes most competitive with pressure increasing until a more stable structure γ-B₂₈ (28 atoms in the unit cell) emerges at 19 GPa.^[14] The γ-B₂₈ structure contains B₁₂ icosahedron and B₂ pairs, which can be considered as “anions” and “cations”, respectively.^[14] A denser α-Ga phase of poor metallicity appears at 89 GPa and undergoes pressure-induced metal–insulator transitions above 300 GPa. Moreover, this phase is predicted to be superconductor with T_c value of 6 K at 175 GPa.^[15] Phosphorus, which has five valence electrons per atom, shows a fascinating structural phase transition under pressure and exhibits some unique phases. Under ambient conditions, black P with Cmca phase appears as a semiconductor structure and transforms into a semimetallic phase (R-3m) at 4.5 GPa. The R-3m phase undergoes a structural transition to a metallic SC phase (space group Pm3m) at 10 GPa and it is stable until 107 GPa.^{[16, 17]} Recently, superconductivity in the SC phase has been widely explored, and T_c value reaches 4.5–13 K under high pressure; however, the dependence of T_c on pressure remains controversial.^{[18, 19]} Then P transforms into the simple hexagonal phase of P6/mmm at 132 GPa via an incommensurate modulation Cmmm,^[20–22] and at higher pressure, another Im-3m phase becomes the most stable structure.^[22] These interesting properties of B and P have aroused strong curiosity about exploring boron–phosphorous under high pressure.

In recent decades, BP and B₆P have been investigated theoretically and experimentally. BP with F-43m structure can be synthesized experimentally from ambient environment to 110 GPa.^[23] A recent theoretical study predicted that the zinc-blende structure transforms into rocksalt structure at 142 GPa.^[24] The other theoretical result shows that the zincblende transforms to the nickel arsenide at 133.26 GPa, and then transforms to rocksalt at 211.99 GPa.^[25] However, Zhang et al. reported that F-43m structure transforms to novel C2/c structure at about 113 GPa and then changes to P4₂/mnm above 208 GPa.^[26] As for B₆P, a recent experiment on Raman about B₆P shows that the R-3m structure transforms to a distorted structure at 80 GPa.^[27] In this work, we use the evolutionary local random structural prediction method along with first-principles calculations to investigate structures and properties of B–P compounds in the pressure range of 0–200 GPa. Only BP and B₆P compounds are thermodynamically stable under high pressures. Electronic properties, electronic band structures, bonding patterns, and dynamical stability of F-43m BP and R-3m B₆P are systematically investigated. Based on the experimental x-ray diffractions and theoretical simulations, we have proved the correctness of the structures we proposed. The infrared (IR) and Raman spectroscopies will provide reliable data for future scientific research about B–P compounds. Meanwhile, the decomposing mechanism of BP and B₆P is deeply analyzed.

The structure prediction for B–P compounds at high pressures is performed by the in-house developed evolutional local random (ELocR) code.^[28–30] The fully structural relaxation, electronic band structures, electronic localization functions (ELF), charge density difference, and Bader charge analysis are carried out by the Vienna ab initio simulation package (VASP).^[31] We choose the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) for the exchange–correlation function. The projector augmented wave (PAW) method is selected, s²p¹ and s²p³ are considered as valence electrons for B and P, respectively. The plane-wave cutoff energy of 950 eV is adopted to ensure the convergence of total energy calculations. The electronic self-consistent calculation is stopped until the energy convergence is less than 10⁻⁶ eV.

Phonon calculations are carried out by a finite displacement approach through the Phonopy code.^[32] The zero-point energy (ZPE) is obtained by

The Raman and IR spectra are calculated by Cambridge serial total energy package (CASTEP) code.^[33] We use the PBE and GGA for the exchange–correlation function with norm-conserving pseudopotentials in the code, and the Brillouin zone sampling grid with a spacing of 2π ×0.03 Å⁻¹.

To ensure the accuracy of our calculation, we have tested the convergence of energy cut-off and grid density, as shown in Fig. S1 in the Supplementary Material. The energy reaches convergence when the cut-off energy and grid density are 950 eV and 2π ×0.03 Å⁻¹ (i.e., the grid size of 16×16×16 in the BP at 50 GPa in Fig. S1), respectively.

Variable-composition structure searches for a variety of B–P stoichiometries are explored in the pressure range of 0–200 GPa. Whether these structures of different stoichiometries are thermodynamically stable depends on their formation enthalpies. The formation enthalpies is determined by the equation

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Calculated formation enthalpies of B–P system with respect to B and P at selected pressures. Solid symbols (solid lines) represent the thermodynamically stable stoichiometries at corresponding pressures, and the open symbols (dashed lines) represent metastable ones.

The enthalpy–pressure curves of BP and B₆P are calculated in the range of 0–200 GPa. According to the enthalpy curves in Fig. 2(a), the cubic F-43m phase (Fig. 3(a)) is the most stable structure below 85 GPa. And the predicted phase is in good agreement with previous theoretical and experimental results.^[34–36] In a recent theoretical study, C2/m BP is proposed a stable metallic phase with a superconducting critical temperature of 9.4 K–11.5 K at 113 GPa.^[26] Unfortunately, C2/m BP would be decomposed into B and P above 85 GPa based on our calculation of formation enthalpies. The enthalpy–pressure curves of several possible candidate structures of B₆P have also been explored. As shown in Fig. 2(b), the trigonal structure of R-3m (Fig. 3(b)) which is proposed by Guo et al,^[37] is stable below 81.4 GPa. The boron icosahedron is discovered in the R-3m structure and connects with P atoms. In the icosahedron, three B atoms are located in the upper part, three B atoms are in the lower part, and the other six equivalent B atoms form the middle part. The structural information for F-43m and R-3m phases are summarized in Table S1. Using the radiational wavelengths of 1.542 Å and 1.789 Å, the x-ray diffraction (XRD) patterns of the predicted BP and B₆P are simulated and compared with powder diffraction file (PDF) standard card at ambient pressure. As shown in Fig. S2, the calculated XRD patterns agree well with the PDF standard card.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Calculated enthalpies per atom as a function of pressure. The formation enthalpies of (a) BP and (b) B6P are shown relative to the decomposition to B+P and 6B+P, respectively.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Crystal structures of (a) F-43m BP and (b) R-3m B6P at 40 GPa. Purple and green spheres depict P and B atoms, respectively.

The enthalpy–pressure curves of BP and B₆P are calculated in the range of 0–200 GPa. According to the enthalpy curves in Fig. 2(a), the cubic F-43m phase (Fig. 3(a)) is the most stable structure below 85 GPa. And the predicted phase is in good agreement with previous theoretical and experimental results.^[34–36] In a recent theoretical study, C2/m BP is proposed a stable metallic phase with a superconducting critical temperature of 9.4 K–11.5 K at 113 GPa.^[26] Unfortunately, C2/m BP would be decomposed into B and P above 85 GPa based on our calculation of formation enthalpies. The enthalpy–pressure curves of several possible candidate structures of B₆P have also been explored. As shown in Fig. 2(b), the trigonal structure of R-3m (Fig. 3(b)) which is proposed by Guo et al,^[37] is stable below 81.4 GPa. The boron icosahedron is discovered in the R-3m structure and connects with P atoms. In the icosahedron, three B atoms are located in the upper part, three B atoms are in the lower part, and the other six equivalent B atoms form the middle part. The structural information for F-43m and R-3m phases are summarized in Table S1. Using the radiational wavelengths of 1.542 Å and 1.789 Å, the x-ray diffraction (XRD) patterns of the predicted BP and B₆P are simulated and compared with powder diffraction file (PDF) standard card at ambient pressure. As shown in Fig. S2, the calculated XRD patterns agree well with the PDF standard card.

The properties of decomposing B–P compounds were then explored by analyzing the formation enthalpies. The different enthalpy curves of B₆P, BP+5B, and 6B+P at 0–200 GPa are displayed in Fig. 4. B₆P is thermodynamically unstable and decomposes into BP and B mixtures above 81.4 GPa, and BP completely becomes B and P mixtures at 85 GPa. Generally, BP and B₆P are thermodynamically stable in lower pressure range, and decompose into B and P at higher pressures. It is necessary to consider the ZPE in B–P compounds due to the lighter elements of B.^[38] Considering the ZPEs of BP, B₆P, B, and P with the harmonic approximation, the shift of decomposition pressure to lower pressures can be found in Fig. S3, i.e., from 85 GPa to 82.8 GPa for BP and 81.4 GPa to 80.1 GPa for B₆P, respectively.

The properties of decomposing B–P compounds were then explored by analyzing the formation enthalpies. The different enthalpy curves of B₆P, BP+5B, and 6B+P at 0–200 GPa are displayed in Fig. 4. B₆P is thermodynamically unstable and decomposes into BP and B mixtures above 81.4 GPa, and BP completely becomes B and P mixtures at 85 GPa. Generally, BP and B₆P are thermodynamically stable in lower pressure range, and decompose into B and P at higher pressures. It is necessary to consider the ZPE in B–P compounds due to the lighter elements of B.^[38] Considering the ZPEs of BP, B₆P, B, and P with the harmonic approximation, the shift of decomposition pressure to lower pressures can be found in Fig. S3, i.e., from 85 GPa to 82.8 GPa for BP and 81.4 GPa to 80.1 GPa for B₆P, respectively.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The formation enthalpies of B6P and BP+5B relative to the decomposition to 6B+P in the pressure range of 0–200 GPa. The inset exhibits in the 50–100 GPa.

We then analyze the structural decomposing of BP and B₆P under high pressures. Enthalpy has two parts: the product of pressure and volume PV, and the internal energy U. The following formulas are adopted to quantify the difference of V and U of BP and B₆P relative to the decomposition products B and P:

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The calculated product of pressure and volume PV and the internal energy U of BP and B6P with respect to B and P as a function of pressure. (a) The and of BP; (b) the and of B6P; (c) and (d) the calculated volume per BP and B6P with respect to the decomposition to B+P and 6B+P at corresponding pressures.

Table 1.

Table 1.

Table 1. The calculated contributions of and to thermodynamic stability of BP and B6P at different pressures. . 0 GPa 20 GPa 40 GPa 60 GPa 80 GPa 100 GPa BP ( ) 100.0% 80.0% 67.5% 58.5% 51.7% 46.3% BP ( ) 0.0% 20.0% 32.5% 41.5% 48.3% 53.7% B6P ( ) 100.0% 79.6% 66.9% 57.9% 51.1% 46.1% B6P ( ) 0% 20.4% 33.1% 42.1% 48.9% 53.9%

Table 1. The calculated contributions of and to thermodynamic stability of BP and B6P at different pressures. .

Meeting the mechanical stability is a fundamental condition for a stable crystal structure, and the elastic constant matrix should be positive depending on Boron–Huang mechanical stability.^[39] The calculated elastic constant matrix C_ij of BP and B₆P at selected pressures are listed in Table S2. For cubic structure of F-43m BP, it satisfies these inequalities: , , |C₁₂|, and , indicating the mechanical stable property.^[40] For R-3m B₆P, it also satisfies the inequalities: , |C₁₂|, and ( . Therefore, the F-43m BP and R-3m B₆P are mechanically stable at corresponding pressures.^[40]

We have also evaluated the hardness of BP and B₆P with the equation

The hardness of B₆P is 42.3 GPa at ambient pressure, and it is a potential superhard material with H_v higher than 40 GPa. For BP, which is the same group boride as the c-BN, it is also crystallized in cubic crystal; however, the H_v of BP is 36.1 GPa, which is smaller than the value of 40.9 GPa in c-BN. The calculated Pugh's ratios of BP and c-BN are 0.99 and 0.79, and the G values are 162.28 for BP and 315 for c-BN,^[42] respectively. The difference of H_v is attributed to the shear modulus G.

To determine the dynamical stability of the crystal structure, we have calculated the phonon band structure and projected phonon density of states (PHDOS) of BP and B₆P at selected pressures, as shown in Fig. 6. At 0 GPa and 40 GPa, both structures demonstrate their dynamical stability without negative frequency in the entire Brillouin zone. From Figs. 6(a) and 6(b), we can see that the low frequency modes mainly come from the vibrations of P atoms, and the higher frequency modes are mainly associated with B atoms due to the heavier atomic mass of P than that of B. Due to the higher content of atoms of B than that of P in the unit cell of R-3m B₆P, the contribution to vibration modes mainly comes from B atoms (Figs. 6(c) and 6(d)). With the pressure increasing, the phonon modes tend to harden in both F-43m BP and R-3m B₆P.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Phonon dispersions and projected phonon density of states for (a) F-43m BP at 0 GPa, (b) F-43m BP at 40 GPa, (c) R-3m B6P at 0 GPa, and (d) R-3m B6P at 40 GPa.

Raman spectra, which provide important information in the research of crystal structures and bonding properties, are simulated and compared with experimental results for BP and B₆P. For R-3m B₆P phase, the Raman active frequency of our theoretical calculation is in good agreement with previous experimental result at 0 GPa.^[27] With the pressure increasing, B₆P shows an expected phenomenon of a hardening of almost all the Raman modes, and a Raman mode ascending to a very high frequency of 1438.83 cm⁻¹ at 80 GPa. There is no Raman mode below 200 cm^−1[27] till the decomposition of B₆P compound. For the F-43m BP, there is only a Raman mode of 781.30 cm⁻¹ at atmospheric pressure, which is slightly lower than 794 cm⁻¹ and 823 cm⁻¹ by Vladimir et al.^[16] The relationship between Raman mode and pressure is depicted in Fig. S4. The Raman peak of F-43m phase tends to a higher frequency due to the shrinkage of B–P bond as the pressure increases. Our calculation is consistent with Vladimir's experimental result on the evolution of Raman shift with pressures.^[16] Furthermore, the vibrational modes of IR and Raman on the phases of F-43m BP and R-3m B₆P can be classified by the irreducible representation of the point groups T_d and D_3d, as shown in Table S3. There are two atoms in the primitive cell of the F-43m structure; therefore, six vibrational modes are produced including three optical modes 3T₂ which are third-degree degeneration. So , and T₂ mode is both IR active and Raman active. In R-3m phase, there are 14 atoms in the primitive cell, and therefore, there are 42 vibrational modes including 39 optical modes, arrived at . The and E_u modes are IR active; E_g and A_1g modes are Raman active; A_1u and A_2g are silent; E_u and E_g are second-degree degeneration.

In the meantime, the electronic properties of BP and B₆P are calculated. The band structures of BP and B₆P near Fermi energy at high pressure are shown in Fig. S5. Obviously, BP and B₆P are indirect bandgap semiconductors, with the feature that the maximum value of the valence band and the minimum value of conduction band are not on the same high symmetry point. The electron localization function and charge density difference are calculated to investigate the bonding characteristics and distribution of electrons in B–P compounds. The ELF of 1 represents electron complete localization, the ELF of 0.5 reflects the probable electron-gas-like pairs, and the value of 0 expresses non-electron localization. Figure S6(a) and S6(b) disclose three dimensional ELFs for F-43m BP and R-3m B₆P at 40 GPa with an isosurface value of 0.75, respectively. Apparently, a large number of valence electrons accumulate in the interstitial regions of both crystals. According to Figs. S6(c) and S6(d), the valence electrons are localized between B and P atoms, indicating that B and P have strong covalent interaction for F-43m phase. Meanwhile, valence electrons are found to localize between the neighbouring B atoms and P atoms in R-3m. In R-3m phase, there are two kinds of B–B bonds, one is located on the boron icosahedra surface (B–B¹) and the other connects the icosahedra (B–B²). The ELF value of B–B¹ bond reaches 0.84, which is much smaller than the value of 0.97 of the B–B² bond. In addition, the ELF values of the B–P and P–P bond are 0.96 and 0.93, respectively. The value of ELF around the outside of boron icosahedra is much larger than that on the surface of boron icosahedra. The relationship between the ELF values and pressure for all bonds of B–P compounds is shown in Table S4. Additionally, in the F-43m structure at 0 GPa, Bader charge analysis discloses the charge transfer from B to P. Due to the variety of bonds in R-3m B₆P, the result based on Bader charge analysis becomes complex. Although the total charge transfer from B to P can be found in R-3m phase, the transfers of and are also observed here, as shown in Table S5.

In summary, we have executed systematic searching for crystal structures of B–P system in the pressure range of 0–200 GPa by utilizing evolutionary structure searches. Two stoichiometries of BP and B₆P are discovered to be thermodynamically stable below 85 GPa and 81.4 GPa, respectively. The electronic energy band structures demonstrate that both BP and B₆P are indirect bandgap semiconductors at corresponding pressure interval. The analyses of ELF display the complicated bonding characteristic of BP and B₆P. With further investigation of IR and Raman spectra of BP and B₆P, we have observed that BP and B₆P express a phenomenon with a hardening of all the Raman modes. Accurate vibrational modes are helpful in further experimental study about interesting phenomena of B–P compounds at high pressures. The analyses of decomposing mechanism of BP and B₆P show that the decomposing of BP and B₆P is dominated by the volumes of B, P, and corresponding B–P compounds. Less-packed BP and B₆P would be more likely to be thermodynamically unstable under high pressure. In addition, our conclusions will provide a great push to research to explore the structures and properties of other borides.

Parts of the calculations were performed in the High Performance Computing Center (HPCC) of Jilin University.